Statistical Physics for Humanities: 
A Tutorial 

Dietrich Stauffcr 

The image of physics is connected with simple "mechanical" deterministic events: that 
an apple always falls down, that force equals mass times acceleration. Indeed, ap- 
plications of such concept to social or historical problems go back two centuries (pop- 
ulation growth and stabilisation, by Malthus and by Verhulst) and use "differential 
equations", as recently revierwed by Vitanov and Ausloos [2011]. However, since even 
today's computers cannot follow the motion of all air molecules within one cubic cen- 
timeter, the probabilistic approach has become fashionable since Ludwig Boltzmann 
invented Statistical Physics in the 19th century. Computer simulations in Statistical 
Physics deal with single particles, a method called agent-based modelling in fields 
which adopted it later. Particularly simple are binary models where each particle has 
only two choices, called spin up and spin down by physicists, bit zero and bit one by 
computer scientists, and voters for the Republicans or for the Democrats in American 
politics (where one human is simulated as one particle). Neighbouring particles may 
influence each other, and the Ising model of 1925 is the best-studied example of such 
models. This text will explain to the reader how to program the Ising model on a 
square lattice (in Fortran language); starting from there the readers can build their 
own computer programs. Some applications of Statistical Physics outside the natural 
sciences will be listed. 

1 Introduction 

Learning by Doing is the intention of this tutorial: readers should learn how 
to construct their own models and to program them, not learn about the great 
works of the author [Stauffer et al 2006] and the lesser works of his competitors 
[Billari et al 2006]. 

Already Empedokles is reported to have 25 centuries ago compared humans to 
fluids: Some are easy to mix, like wine and water; and some, like oil and water, 
refuse to mix. Newspapers can give you recent human examples, and physicists 
and others have taken up the challenge to study selected problems of history 
and other humanities [Castellano et al 2009] with methods similar to physics. 
In the opposite direction, Prados [2009] uses the physics dream of "unified field 
theory" to describe his history of the Vietnam war. 

The next section will recommend ways how to construct models, the fol- 
lowing one how to program a simple Ising model on a square lattice, and a 
concluding section will list some applications. An appendix will introduce the 
Fortran programming language. 
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2 Model Building 



2.1 What is a Model? 

"Models" in physics and in this tutorial usually deal with the single elements of 
a system and how their interactions produce the behaviour of the whole system. 
Outside of physics, a model may just be any mathematical law approximating 
reality. Thus the statement that human adult mortality increases exponentially 
with age is often called the Gompertz model by demographers but the Gompertz 
rule by physicists; the latter ones use the Penna model of individuals undergo- 
ing genetic mutations and Darwinian selection to simulate a large population 
perhaps obeying Gompertz [Stauffer et al 2006]. 

2.2 Binary versus more complicated models 

If no previous work on a general model is known, I recommend to start with 
binary variables where each particle has only two choices, called spin up and 
spin down by physicists, bit zero and bit one by computer scientists, occupied 
or empty for percolation, and voters for the Republicans or for the Democrats 
in American politics (where one human is simulated as one particle) . Of course 
reality is more complicated, but we want to understand reality: Does it agree 
with the simplest possible model? (Simulations for pilot training etc. are differ- 
ent [Bridson and Batty 2010]). Thus we follow the opinion of Albert Einstein 
that a model should be as simple as possible, but not simpler. If you want to 
simulate traffic jams [Chowdhury et al. 2000] in cities, the colour of the cars is 
quite irrelevant, but for visibility in the dark the colour matters. 

Once the binary case has been studied, one can go to more than two choices. 
If three groups are fighting each other, obviously three choices are needed in 
a model [Lim 2007]. Even in a two-party political system like the USA, other 
candidates were important in Florida for the US presidential elections of 2000. 
The opinions of people [Malarz et al 2011] are in reality continuous and can be 
modelled by one or several real numbers between zero and unity, or between 
minus infinity and plus infinity. Nevertheless it is standard practice in opinion 
polls to allow only a few choices like full agreement, partial agreement, neutral- 
ity, partial disagreement, full disagreement. And in elections one can only vote 
among the discrete number of candidates or parties which are on the ballot. In 
the Kosovo opinion of the International Court of Justice (July 22, 2010) one 
judge criticised the binary tradition of either legal or illegal, stating that toler- 
able is in between; nevertheless the court majority stayed with the binary logic 
of not illegal. 

Physicists like to call the binary variables "spins" but readers from outside 
physics should refrain from studying the quite complicated spin concept in quan- 
tum mechanics. Spins are simply up and down (1 or 0; 1 or -1). Similary, don't 
be deterred if physicist talk about a Hamiltonian; in most cases this is just the 
energy from high school. 
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2.3 Humans are neither Spins nor Atoms 

Of course, that is true, but it docs not exclude that humans arc modelled like 
spins. Modern medicine blurred the boundary between life and death; never- 
theless we usually talk about people having been born in a certain year, and 
having died in another year, as if they were binary up-down variables. Reasons 
for death are complicated and I don't even know mine yet; nevertheless de- 
mographers like Gompertz estimated probabilities for dying at some age. Such 
probabilities are rather useless for predicting the death of one individual, but 
averaged over many people they may give quite accurate results. When I throw 
one coin I do not know how it will fall; when I throw thousand coins, usually 
about half of them fall on one side and the others on the other side (law of large 
numbers). If I throw 1000 coins and all show "head", most likely I cheated. 
Thus to simulate one person's opinion and decisions on a computer docs not 
seem to be realistic; to do the same for millions of people may give good aver- 
age properties, like the number of deaths at the age of 80 to 81 years, or the 
fraction of voters selecting the parties in an upcoming election. The whole in- 
surance industry is based on this law of large numbers. Humans are not spins 
but many humans together might be studied well by spin models. 

2.4 Deterministic or Statistical ? 

Non-physicists often believe that physics deals with deterministic rules: The 
apple falls down from the tree and not up; force equals mass times acceleration; 
etc. (Or they have heard of quantum-mechanical probability and apply that 
to large systems where such quantum effects should be negligibly small.) In 
this sense the cause of World War I was seen as a consequence of the arms 
race modelled by deterministic differential equations for averages [Richardson 
1935]. And the decay of empires was described [Geiss 2008] as starting at the 
geographical periphery, since the influence from the center decreases towards 
zero with increasing distance, just as the gravitational force between the sun 
and its planets; see also [Diamond 1997, epilog]. 

Because of its historical importance let us look into Richardson's papers of 
1935: Two opposing (groups of) nations change their preparedness x for war 
because of three reasons: 1) the war preparedness of the other side; 2) fatigue 
and expense; 3) dissatisfaction with existing peace treaties. Reasons 1 and 3 
increase and reason 2 decreases the war preparedness; reason 1 is proportional 
to the x of the opponent, reason 2 to the own x, and reason 3 independent of 
x\ and X2- Even complete disarmament [x\ = x^ = 0) at one moment does 
not help if there is dissatisfaction with the existing peace. And if reason 1 is 
stronger than reason 2 for both sides, both x increase exponentially with time 
towards infinity. His second paper details the mathematical solutions for these 
linear coupled inhomogenous differential equations. 

More recently, the widespread use of computers shifted the emphasis to more 
realistic probabilistic models, using random numbers to simulate the throwing 
of coins or other statistical methods. This "Statistical Physics" and a simple 
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example are the subject of the next section. 

3 Statistical Physics and the Ising Model 

3.1 Boltzmann Distribution 

No present computer can simulate the motion of all air molecules in a cubic 
centimeter. Fortunately, Ludwig Boltzmann about 150 years ago invented a 
simple rule. The molecules move at temperature T with a velocity which can 
change all the time but follows a statistical distribution: The probability for a 
velocity v is proportional to exp(—E/T), where E is the kinetic energy of the 
molecule due to its velocity. The same principle applied to a binary choice, 
where a particle can be in two states A and B with energies Ea and Eb means 
that the two probabilities are 

PA = |exp(-£7 A /T); p B = \ cxp(-£ B /T); (la) 

Z = cM-Ea/T) + cxp(-E B /T) (lb) 

since the sum over all probabilities must be unity. More generally, a configura- 
tion with energy E is in thermal equilibrium found with probability 

p=^exp(-E/T); Z = ]T exp(-£7/T) (2) 

where the sum runs over all possible states of the system. Z is called the 
partition function, one of the rare cases where the German word for it, Zu- 
standssumme = sum over states, is clearer and shorter. The temperature T 
is measured neither in Celsius (centigrade) nor Fahrenheit but T — at the 
absolute zero temperature (about -273 Celsius below the freezing temperature 
of water) and moreover is measured in energy units. (If T is measured in Kelvin, 
the corresponding energy is k B T where k B is the Boltzmann constant and set to 
unity in the present tutorial.) The function exp(x) is the exponential function, 
also written as e x , which for integer x means the product of x factors e ~ 
2.71828; e x = 2 X / - 69315 = io°- 4343a: . 

Eqs.(l,2) can be regarded as axioms on which Statistical Physics is built, 
like the Parallel Axiom of Euclidean Geometry, but in some cases they can be 
derived from other principles. Humanities are allowed to use these ideas of 
Boltzmann since history institutes were named after him. 

3.2 Ising Model 

In 1925, Ernst Ising (born in the heart of the city this author lives in) finished 
his doctoral dissertation on a model for ferromagnetism, which became famous 
two decades later and was shown to apply to liquid- vapour transitions half a 
century later. We assume that each site of a lattice (e.g. a square lattice where 
each site i has four neighbours: clockwise up, right, down, left) carries a spin 
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Si = ±1 (up or down). Neighbouring spins "want" to be parallel, i.e. they have 
an energy — J if they are in the same state and an energy + J if they are in the 
two different states. Moreover, a "magnetic" field H between minus infinity and 
plus infinity (also called B) tries to orient the spins in its own direction. The 
total energy then is 

E = -Jj^ SiSj-Hj^Si (3) 

<ij> i 

where the first sum goes over all ordered pairs of neighbor sites i and j. Thus 
the "bond" between sites A and B appears only once in this sum, and not twice 
(for i = A, j = B as well as for i = B, j = A). The second sum runs over all 
sites of the system. Thus 2 J is the energy to break one bond, and 2H is the 
energy to flip a spin from the direction of the field into the opposite direction. 
As discussed before in the Boltzmann subsection, the higher the energy E is 
the lower is the probability to observe this spin configuration; at infinitely high 
temperatures T all configurations are equally probable; at T = all spins must 
be parallel to each other and to the field H in equilibrium. The "magnetisation" 
M is the number of up spins minus the number of down spins, 

M = Y,Si . (4) 

i 

Computer simulations of this Ising model will be described in the appendix. 

Applied to human beings, this Ising model could represent two possible 
opinions in a population; everybody tries to convince the neighbours of the 
own opinion (J), and in addition the government H may try to convince the 
whole population of its own opinion. The temperature then gives the tendency 
of the individuals not to think like the majority of their neighbours and the 
government. Zero temperature thus means complete conformity, and infinite 
temperature completely random opinions. 

Theories with paper and pencil in two dimensions as well as computer simu- 
lations give M(T, H). In particular, for H — and in more than one dimension, 
the equilibrium magnetisation M = ±Mq is a non-zero spontaneous magneti- 
sation for T < T c and is zero for T > T c where T c is the critical or Curie 
temperature (named after Pierre Curie, not his more famous wife Marie Curie). 
On the square lattice, T c /J ~ 2.27 is known exactly, in three dimensions only 
numerically; in one dimension there is no transition to a spontaneous magneti- 
sation, as Ernst Ising had shown: T c = 0. 

The above model obeys Isaac Newton's law actio = —reactio: The sun 
attracts the earth with the same force as the earth attracts the sun, only in 
opposite direction. Human relations can also be unsymmetrical: He loves her 
but she does not love him. Then the bond between sites i and j may be directed 
instead of the usual case of an undirected bond. For example, if i influences j 
but j does not influence i, flipping the spin at j from parallel (Sj = +5,) to 
antiparallel (Sj = —Si) to the spin at i can cost an energy 2J while flipping Si 
at constant Sj costs nothing. In this case no unique energy E(Si, Sj) is defined 
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for this spin pair and thus such models have been much less studied in the 
physics literature. (In this example, one could gain a lot of energy from nothing 
by the cyclic process of flipping j from antiparallel to parallel, gaining energy 
2 J, then flipping i from parallel to antiparallel, costing nothing, then flipping j 
again and so on. Such a perpetuum mobile does not exist in physics.) 

Outside physics, of course, one can commit such crimes against energy con- 
servation and forget all probabilities proportional to exp(—E/T). Instead one 
can assume arbitrary probabilities, as long as their sum equals one. For example, 
if an element has three possible states A, B and C, then one may assume that 
with probability p an A becomes B, a B becomes C, and a C becomes A; with 
probability 1—p the element does not change, and there is no backward process 
from C to B to A to C. Then one has a circular perpetuum mobile, which may 
be realistic for some social processes. Physicists sometimes distinguish between 
dynamics (when the changes are determined fully by energy or force) and ki- 
netics (when additional assumptions, like the probabilites of eqs.(l) are made); 
probabilities independent of energy/force are then kinetics, as is most of the ma- 
terial described here. And usually to find a stationary or static equilibrium, one 
has to wait for many non-equiilibrium iterations in the simulation (in a static 
situation, nothing moves anymore; in a stationary simulation the averages are 
nearly constant since changes in one direction are mostly cancelled by changes 
of other elements in the opposite direction. There are many choices, one needs 
not mathematics, but mathematical thinking: precise and step-by-step. 

3.3 Warning against Mean Field Approximation 

This subsection contains many formulas and can be skipped. If you want to 
get answers by paper and pencil, you can use the mean field approximation 
(also called molecular field approximation), which in economics correponds to 
the approximation by representative agent. Approximate in the first sum of 
Eq.(3) the Sj by its average value, which is just the normalised magnetisation 
to = M/L 2 = £\ Si/L 2 . Then the energy is 

E = - j 5 Sim ~ H ? St = ~ Heff 2 Si 

<ij> i 

with the effective field 

H e ff = H + ^ m = H + qm 

3 

where the latter sum runs over the q neighbours only and is proportional to the 
magnetization to. Thus the energy Ei of spin i no longer is coupled to other 
spins j and equals ±H e ff. The probabilities p for up and down orientations, 
according to Eq.(l), are now 

P(Si = +1) = \ cMHeff/T); p{S t = -!) = ! exp(-if e// /T) 
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and thus 



m = p(Si = +l)-p(S l = -1) = tanh(H e// /T) = tanh[(i? + qm)/T] 

with the function tanh(x) = (e x — e~ x )/(e x + e~ x ). This implicit equation can 
be solved graphically; for small m and H/T, tanh(x) = x — x 3 /3 + . . . gives 

H/T = (1 - T c /T)m + im 3 + . . . ; T c = qJ 

related to Lev Davidovich Landau's theory of 1937 for critical phenomena (T 
near T c , m and H/T small) near phase transitions. 

All this looks very nice except that it is wrong: In the one-dimensional Ising 
model, T c is zero instead of the mean field approximation T c = qj. The larger 
the number of neighbours and the dimensionality of the lattice is, the more ac- 
curate is the mean field approximation. Basically, the approximation to replace 
SiSj by an average Sim takes into account the influenve of Sj on Si but not 
the fact that this Si again influences Sj creating a feedback. If the mathematics 
of this subsection looks deterrent, just ignore it; you are recommended to use 
computer simulations of single interacting spins, and not mean field theories. 
Outside of physics such simulations are often called "agent based" [Billari et al 
2006, Bonabcau 2002]; presumably the first one was the Metropolis algorithm 
published in 1953 by the group of Edward Teller, who is historically known from 
the US hydrogen bomb and Strategic Defense Initiative (Star Wars, SDI). 

4 Applications 

4.1 Schelling Model for Social Segregation 

Economics nobel laureate Thomas C. Schelling advised the US government on 
war and peace in the 1960s and published a methodologically crucial paper 
(cited more than 500 times a decade later) [Schelling 1971, Fossett 2011, Henry 
et al 2011] which introduced methods of statistical physics to sociology. Each 
Ising spin corresponds to one of two ethnic groups (black and white in big US 
cities) both of which prefer not to be surrounded by the other group. The 
above Ising model then shows that for T < T c segregation emerges without 
any outside control: The simulated region becomes mostly white or mostly 
black. Unfortunately, Schelling simulated a more complicated version of the 
Ising lattice at T = which did not give large "ghettos" like Harlem in New 
York City, but only small clusters of predominantly white and black residences. 
Only by additional randomness [Jones 1985] (cited only 8 times, mostly by 
physicists) can "infinitely" large ghettos appear. It is easier to just simulate 
the above Ising model [Sumour et al 2008,2011], which also allows people of 
one group to move into another city (or away from the simulated region) to be 
replaced by residents of the other group. Sumour et al also cite earlier Schelling- 
type simulations by physicists. Fore nearly three decades physicists ignored the 
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Schclling model; now the sociologists ignore the physics simulations of the last 
decade and the much earlier Jones [1985] paper. 

As examples we give two ilsing-model figures from Miiller et al [2008] where 
people also increase their amount T of tolerance = social temperature if they see 
that their whole neighbourhood belongs to the same group as they themselves. 
And afterwards they slowly forget and reduce their tolerance. With changing 
frogeting rate one may either observe lots of small clusters, Fig.l, or one big 
ghetto, Fig. 2. 



Small clusters in modified Ising model 




50 100 150 200 250 300 350 400 



Figure 1: Cluster formation for slow forgetting in the Ising modification of 
Miiller ct al [2008] ; no large ghetto is formed. 

About simultaneously with this Schelling paper, physicist Weidlich started 
his sociodynamics approach to apply the style of physics to social questions 
[Weidlich 2000], see also [Galam 2008]. 

4.2 Sociophysics and Networks 

A good overall review of the sociophysics field (of which the Schelling model is 
just one example) was given by Castellano et al [2009]. Of particular interest 
is the reproduction of universal properties of election results with many candi- 
dates: The curves of how many candidates got n votes each are similar to each 
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Modified Ising model: feedback, memory loss 
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Figure 2: In the same model of the previous figure, a large ghetto is formed if 
people forget faster their learned tolerance. 

other [Castellano 2009]. Car traffice [Chowdhury et al [2000], economic mar- 
kets [Bouchaud and Potters 2010, Bonabeau 2002], opinion dynamics [Malarz et 
al 2011], stone-age culture [Shennan 2001], histophysics [Lam 2002], languages 
[Schulze et al 2008] , Napoleon's decision before the battle of Waterloo [Mongin 
2008], religion [Ausloos and Petroni 2009], political secession from a state [Lu- 
stick 2011], demography on social networks [Fent et al. 2011], insurgent wars 
in Iraq and Afghanistan [Johnson et al. 1944], ... are other applications. For 
example, non-physicists [Holman et al 2011] have acknowledged that physicists 
Serva and Petroni (2008) may have been the first to use Levinstein distances to 
calculate the ages of language groups. (These distances are differences between 
words for the same meaning in different languages.) 

Students in class may sit on a square lattice, but normally humans do not. 
(In universities, mostly only a part of the lattice is occupied, which physicists 
simulate as a "dilute" square lattice.) They may be connected by friendship or 
job not only with nearest neighbors but also with people further away. Such 
networks, investigated for a long time by sociologists [Stegbauer and Haeussling 
2010], were studied by physicists intensively for a dozen years [Barabasi 2002, 
Albert and Barabasi 2002, Bornholdt and Schuster 2003, Cohen and Havlin 
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2010]. In the Watts-Strogatz (or "small world") network a random fraction of 
nearest-neighbour bonds is replaced by bonds with sites further away, selected 
randomly from the whole lattice. If that fraction approaches unity one obtains 
the Erdos-Renyi networks, a limit of percolation theory [Flory 1941]. More 
realistic are the scale-free Barabasi- Albert networks, where the network starts 
with a small core and then each newly added site forms a bond with a randomly 
selected already existing member of the network. The selection probability is 
proportional to the number of bonds which the old member already has acquired: 
Famous people get more attention and more "friends" than others; no lattice 
is assumed here anymore. In all these networks, the average number of bonds 
needed to connect two randomly selected sites increases logarithmically with the 
number of sites in the network, whereas for d-dimensional lattices this average 
number of bonds increases stronger with a power law, exponent 1/d. Having 
simulated one network, one can also study connected sets of networks or other 
social networks [Watts, Dodds, Newman 2002], or demography on them [Fent 
etal2011]. 

The latest application is Statistical Justice: In May 2011, John Demjanjuk 
was sentenced for having helped in 1943 in the murder of more than 28,000 
Dutch Jews in the Nazi concentration camp of Sobibor. One knows who was 
deported but not who survived the transport from the Netherlands to Poland. 
And one does not know which duties the accused had there on which day. Thus 
all acts of the camp guards were regarded as having helped in their murder, and 
the number near 28,000 was estimated from the average death rate during the 
transports. The verdict thus gave neither the name of a murder victim nor the 
day of a murder, but was entirely based on statistical averages. [Times 2011]. 

5 Appendix: How to Program the Ising Model 

The following Fortran manual and program are both short and should encourage 
the reader to learn this technique. 

5.1 Fortran Manual 

Fortran = formula translator an early language (above machine code or as- 
sembler) for computer programming; many others followed and in particular 
C ++ is widespread, but nevertheless this tutorial uses Fortran which is closer 
to plain English and allows to easily find a typical programming error (using an 
array outside its dchncd bounds). If your Fortran program is called name.f, it 
can be compiled with f95 -0 name.f (or 77 instead of 95; = optimisation), 
and executed with ./a. out (or just a. out). The just mentioned error message 
appears when using f95 -f bounds-check name.f ; ./a. out, but execution 
then is much slower and thus instead -0 should be used after error correction. 

Fortran commands usually start in column 7 and end before column 73. A C 
in column 1 signifies a comment for the reader, to be ignored by the computer. 
In column 6 we write a 1 if this line is a continuation of the previous line, 
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while columns 2 to 5 are reserved for labels, i.e. numbers to control the flow of 
commands. For example, GOTO 7 means to jump to the line labelled by 7. 

Variable names start with a letter; names starting with I, J, K, L, M, N 
signify integers without rounding errors, other names are real (floating-point) 
numbers and nearly always have rounding errors. The operations +,—,*,/ and 
SQRT, COS, SIN, EXP etc have their usual meaning except that N/M is always 
rounded downwards to an integer value; e.g. 3/5 is zero. Also I = X means 
rounding downwards. Since the natural logarithm is not an integer it is denoted 
by ALOG instead of log. 

Decisions are made automatically, e.g. 
IF(A.GT.O) B = SQRT (A) 
where .GT. means greater than, with analogous meanings for . LT . , . GE , . LE . , 
EQ . , . NE . , .NOT. , .AND., .OR. 

A loop is executed by 
DO 99 K=M,N 

which means that all lines from this line down to and including the line with 
label 99 are executed for k = m, m + 1, m + 2, n. One may put inner loops 
into outer loops, if needed. 

Arrays need to be declared at the beginning of a program, for example 
through 

DIMENSION A(100, 100) , B(100) , C(L) 

Here, if C has an arbitrary dimension L, then L must be given a value before 
this dimension statement through 

PARAMETER (L = 100) 
and must not be changed throughout the program. Similarly, variables can be 
initialised via a data line like 

DATA L/100/, B/100*1.0/ 
but only once at the beginning of the program, not later again. 
Results are best printed out through 

PRINT *, x, y, z 
Thereafter execution should stop with a STOP line, followed by an END line. 
The statement 

n = n+1 

is not an equality (which then could be simplified to the nonsensical = 1) but a 
command to the computer: to find the place in the memory where the variable n 
is stored, to get the value of n from there, to add one to it, and to store the sum 
in that same memory place as the new value for n. Some computer languages 
therefore use := instead of the simpler but misleading = sign. 

Normally it does not matter whether or not CAPITAL letters are used. The 
computer language Basic is rather similar to Fortran. Now we bring a complete 
program to simulate the Ising model on the square lattice. 

5.2 Ising Model Program 

c heat bath 2D Ising in a field 
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parameter (L=1001,Lmax=(L+2)*L) 
dimension is(Lmax) ,ex(-4:4) 

data t , mcstep , iseed/0 . 90 , 1000 , 1/ , h/+0 . 50/ , ex/9*0 . 0/ 

print *, '#', L, mcstep, iseed,t,h 

x=rand(iseed) 

Lpl=L+l 

LspL=L*L+L 

L2pl=2*L+l 

do 1 i=l,Lmax 

1 is(i)=l 

do 2 ie=-4,4,2 

x=exp(-ie*2 . 0*0 . 4406868/t-h) 

2 ex(ie)=x/(1.0+x) 
do 3 mc=l, mcstep 

mag=0 

do 4 i=Lpl,LspL 
c if (i.ne.L2pl) goto 6 

c do 5 j=l,L 

c5 is(j+LspL)=is(j+L) 
6 ie=is(i-l)+is(i+l)+is(i-L)+is(i+L) 
is(i)=l 

if (randO .lt.ex(ie)) is(i)=-l 
4 mag=mag+ i s ( i ) 

c do 7 i=l,L 

c7 is(i)=is(i+L*L) 

3 if (mc.eq. (mc/100) *100) print *, mc, mag 
stop 

end 

The parameter line fixes the size of the L x L square; the sites in it are 
numbered by one index, typewriter style. Thus the right neighbor of site 2L is 
2L + 1 and sits on the left end of the next line: Helical boundary conditions. 
The lower neighbour of site 2L is 2L + L, the upper neighbour is 2L — L, and 
the left neighbor is 2L — 1. The spins in the top and bottom buffer lines (1 . . . L 
and L 2 + L + 1 . . . L 2 + 2L) stay in their initial up orientation; if instead one 
wants periodic boundary conditions in the vertical direction (better to reduce 
boundary influence) , one has to omit the five comment symbols C at and before 
loops 5 and 7. 

The temperature enters the data line in units of T Cl i.e. T = 0.9T C in this 
example; the known value J/T c = 0.44068 ... is used nine lines later. In the 
same line the field, in units of fc^T, is given as 0.5. In the line after the first 
print statement, rand(iseed) initialises the random number generator (see next 
subsection for warning and improvement); a different seed integer gives different 
random numbers. Later randO produces the next "random" number between 
and 1 from the last one in a reproducible but hardly predictable way. Loop 2 
determines the Boltzmann probabilities iex needed for Eqs.(l) and finishes the 
initialisation. Loop 3 makes mcstep iterations (time steps = Monte Carlo steps 
per spin). Loop 4 runs over all spins except those in the two buffer lines and 
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determines the local interaction energy ie as the sum over the four neighbour 
spins. Then we set the spin to +1, and if the conditions of Eqs(l) so require, 
instead it is set to —1. (Generally, a command is executed with a probability p if 
the random number randO is smaller than p.) We print out the magnetisation 
only every hundred time steps in order to avoid too much data on the computer 
screen. The results are plotted in Fig. 3 and show that the temperature is very 
low: Even though it is only 10 percent below the critical temperature T c , the 
magnctisatuion barely changes from its initial value 1002001 and is after 100 
iterations already in equilibrium, also due to the applied field. In three instead 
of two dimensions, without a field, and at temperatures closer to T c longer times 
are needed for equilibration. 
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Figure 3: Top: Results from the listed Ising program. Bottom: Revolutions 

Perhaps you find it more interesting to simulate revolutions, as in Fig. lb, 
where the field h was equal to the fraction of overturned spins (suggestion of 
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Sorin Solomon for Bornholdt-type model). So, what is difficult about computa- 
tions? Do you know a shorter Fortran manual? 

5.3 Random Numbers 

The above rand produces random numbers in an easily programmed way, but 
often this may be slow and/or bad, or the used algorithm is unknown to the 
user. It is better to program random number generation explicitcly. If you 
multiply by hand two nine-digit integers, you may easily predict the first and 
the last digit of the product, but hardly the digits in the middle, except by 
tediously doing the whole multiplication correctly. Similarly, if ibm is a 32-bit 
odd integer, then the product 
ibm=ibm* 16807 

is again an odd integer, and normally requires 46 bits. (A bit = binary digit 
is a zero or one in a computer.) The computer throws away the leading bits 
and leaves the least significant 32 bits. The first bit gives the sign, thus plus 
times plus gives minus in about half the cases, in contrast to what you learned 
in elementary school. The last of these remaining bits is predictably always set 
to one (odd integers) but the leading (most significant) bits are quite random. 
(Actually, your computer may do something very similar when you call randO .) 
More precisely, they are pseudo-random; in order to search for errors one wants 
to get exactly the same random numbers when one repeats a simulation with 
the same seed. 

These random 32-bit integers ibm between -2147483647 and +2147483647 = 
2 31 — 1 can be transformed into real numbers through ran=f actor*ibm+0 . 5 
where factor = 0.5/2147483647, but it is more efficient to normalize the 
propabilities p (here = ex/ ( 1 . 0+ex) to the full interval of 32-bit integers through 
(2 . 0*p-l . 0) *2147483647, once at the beginning of the simulation. Then the 
next random integers ibm simply have to be compared with this normalized 
probability. In the above Ising program, one then stores the Boltzmann proba- 
bilities as dimension iex(-4:4) at the beginning, calculates them through 

2 iex(ie)=(2.0*ex/(1.0+ex) - 1 . 0) *2147483647 
in the initialisation, and later merely needs 
ibm=ibm* 16807 

if (ibm. It . iex(ie) ) is(i)=-l 

in the above Ising program. 

With 32 bits the pseudo-random integers are repeated after 2 29 such mul- 
tiplications with 16807, which is a rather small number for today's personal 
computers. It is better to use 64 bits via 

integer*8 ibm, iex 
at the beginning of the program in order to get many more different random 
numbers, using e.g. 

2 iex(ie)=2147483647*(4.0*ex/(1.0+ex) - 2 . 0) *2147483647 
for the normalized 64-bit probabilities. Now the quality is much better without 
much loss in speed; unfortunately one now can make much more programming 
errors involving these random numbers. 
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